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Viscous Effects on the Interaction between the Coplanar 
Decretion Disc and the Neutron Star in Be/X-Ray Binaries 
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ABSTRACT 

We study the viscous effects on the interaction between the coplanar Be-star disc and 
the neutron star in Be/X-ray binaries, using a three-dimensional, smoothed particle 
hydrodynamics code. For simplicity, we assume the Be disc to be isothermal at the 
temperature of half the stellar effective temperature. In order to mimic the gas ejection 
process from the Be star, we inject particles with the Keplerian rotation velocity at 
a radius just outside the star. Both Be star and neutron star are treated as point 
masses. We find that the Be-star disc is effectively truncated if the Shakura-Sunyaev 
viscosity parameter ass ^ 1; which confirms the previous semi- analytical result. In 
the truncated disc, the material decreted from the Be star accumulates, so that the 
disc becomes denser more rapidly than if around an isolated Be star. The resonant 
truncation of the Be disc results in a significant reduction of the amount of gas captured 
by the neutron star and a strong dependence of the mass capture rate on the orbital 
phase. We also find that an eccentric mode is excited in the Be disc through direct 
driving due to a one-armed bar potential of the binary. The strength of the mode 
becomes greater in the case of a smaller viscosity. In a high-resolution simulation with 
c^SS = 0-1 J the eccentric mode is found to precess in a prograde sense. The mass 
capture rate by the neutron star modulates as the mode precesses. 

Key words: accretion, accretion discs - binaries: close - hydrodynamics - instabilities 
- stars: emission-line. Be - X-rays: stars. 



1 INTRODUCTION 



(Stella, Wliite fc Rosner 1986| see also [Negueruela et al 



Tlie Be /X-i ' a^i biiiaiies represeiit ilie laigesi subclass of liigli- 
mass X-ray binaries. About two-thirds of the identified sys- 
tems fall into this category. These systems consist of a Be 
star (i.e., a B star with an equatorial disc) and, gener- 
ally, a neutron star. The orbit is wide (several tens of days 
^ Poib ^ several hundred days) and eccentric (0.1 < e < 0.9). 

Most of the Be/X-ray binaries show only transient X- 
ray activity due to transient accretion of the circumstellar 
matter of the Be star, while some show persistent X-ray 



199S). These features imply a complicated interaction be- 
tween the Be-star envelope and the neutron star. 

A Be star has a two-component extended atmosphere, 
a polar region and a cool (~ 10'' K) equatorial disc. The 
polar region consists of a low-density, fast (~ lO'^kms"') 
outflow emitting UV radiation. The wind structure is well 
explained by the so-called line-driven wind model, in which 
the radiative acceleration results from the scatterin g of the 



stellar radiation in an ensemble of spectral lines (Castor 



emissio 1. Each Be/X-ray binary exhibits some or all of the Abbott & Klein 1975; Abbott 198S). On the other hand 



following three types of X-ray activity: 



(i) periodic (Type I) X-ray outbursts, coinciding with pe- 
riastron passage (Lx ~ lO'^^^"'^ ergs~'), 

(ii) giant (Type II) X-ray outbursts (Lx > lO^'' ergs~^), 
which siliow no cleai uibilal muduIuLiun, 



the equatorial disc, which is geometrically thin and nearly 
Keplerian, consists of a high-density plasma from which the 
optical emission lines and the IR excess arise. The radial ve- 
locity of the disc is smaller than a few km s~' , at least within 



(iii) pcraiatcnt low - luminoaity 

(Lx<10='' ergs-i) 



A - ray 



cmi33ion 



10 stellar radii (Hanuschik 1994 ; Hanuschik 200C ; Water! 
Marlborough 1994|). Although there is no widely accepted 



model for discs aro und Be stars, the viscou s decretion disc 
model proposed by Lee, Saio & Osaki (1991) explain s many 



of th e observ ed features an d thus seems promising (Porter 



E-ma 1: okazakiOelsa.hokkai-s-u. ac.jp 



199£ ; see also Okazaki 2001 ) . In this model, the matter sup- 
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'10 km/s wind 



qualitativ ely agrees with the result found by Zamanov et 



B e -Star- 



rs 



near-Keplenan disk 



with < 1 km/s outflow 



^..1 pcri- 
astron 
eccentric 
nrtnt 



Figiiro 1 ■ Schematic view nf a Be/X-ray binary, taken from 
Okazaki fc Neguerucla (2001b). 



plied from the equatorial surface of the star drifts outwards 
because of the viscous effect and forms the disc. The basic 
latious for vi.scous decretion di.scs are the same as tho.se 



equa 

for visc pus accretion discs, except that the sign of M ("mass 



decretion/accretion rate) is opposite. The boundary condi- 
tions for decretion discs, however, are different from those 
for accretion discs. Therefore, the decretion d isc has a struc - 
ture different from that of the accretion disc (Pringle 1991). 



Until quite recently, models for Type I X-ray outbursts 
in Be/X-ray binaries had assumed a large disc around the 
Be star so that the neutron star can accrete ga s when it 
passes through the di sc nea r periastron. However, [Negueru- 



ela & Okazaki (2001) and Okazaki & Negucruela (2001a 



recently performed a semi-analytical study based on the 
viscous decretion disc model for Be stars and showed that 
the Be disc in Be/X-ray binaries is truncated at a radius 
smaller than the periastron distance, as long as ass ^ 
1, where ass is the Shakura-Sunyaev viscosity parameter 
(see Fig. |^ for a schematic view of a Be/X-ray binary). 
The truncation of the disc is due to the resonant torque 
exerted by the neutron star, which removes the angular 
momentum from the disc. The disc radii they obtained 
for seven particular systems (4U 0115-1-63, V 0332-1-53, 
A 0535+262, EXO 2030-^375, 2S 1417-624, GROJ1008-57, 
and 2S 1417—624) are consistent with the X-ray behaviour 
of those syste ms. Moreover, the result is in agreement with 
the result of Reig, Fabregat & Coe (1997) that there is a 
positive correlation between the orbital size and the maxi- 
mum equivalent width of Ha ever observed in a system, a 
measure of the maximum disc size around the Be star in the 
system. 

The truncation of the Be disc in Be/X-ray binaries is 
not surprising. The resonant interaction is important in var- 
ious contexts even in a fly-by encounter with a perturber. 
In fly-by encounters between a disc galaxy and a point mass 
perturber, the energy is always transported from the disc 
to the perturber through the resonant interaction, except 
f or overhead e ncounters where the energy transfer is small 
(Palmer 1983). In distant encounters between a circumstel- 



lar accretion disc and a perturbing mass with r peri /'"disc ^ 2, 
where rpcri and r-disc are the periastron distance and disc ra- 
dius, respectively, the disc material loses energy and angular 
mom entum to the perturber's orbi t through a resonance fea- 
ture dHall, Clarke fc Pringle 1996| ). 

In the case of Be/X-ray binaries, the surface density of 
the Be disc is expected to increase more rapidly than that 
for isolated Be stars, as a consequence of truncation. This 



al. (2001) that the Be discs in Be/X-ray binaries are about 
twice as dense as those around isolated Be stars. The disc 
may finally become optically thick, and b ecome unstable 
to ra diation-driven warping (Pringle 1996; see also Porter 
Multi-wavelength, long-term monitoring observations 



199S 



of V635 Gas, the optical counterpart of 4U 0115-1-63, re- 
vealed that the Be disc in 4U 0115+63 undergoes a quasi- 
cycli c (3 — 5yr) dynamical evolution (Negueruela et al 



2001): after each disc- loss episode, the disc starts reforming. 



grows until it becomes unstable to warping, and after that a 
Type II outburst occurs. Although a direct link between the 
warped disc and the Type II outburst is still missing, the 
dynamical evolution of the Be disc is likely to be the agent 
that controls the X-ray behaviour of the system. 

This way, the truncated disc model, at least qualita- 
tively, explains many of the observed features of Be/X-ray 
binaries. The semi-an alytical model adopted by Negueruela 



Okazaki (2001) and Okazaki & Negueruela (2001a), how- 



ever, only compares the resonant torque integrated over the 
whole orbit with the viscous torque to determine at which 
radius the disc is truncated. Hence, it cannot make a quanti- 
tative prediction about how perfect or imperfect the trunca- 
tion is. Moreover, it predicts nothing about phase-dependent 
features, such as the disc deformation and the change in the 
mass capture rate. 

Therefore, in order to study the efficiency of the res- 
onant truncation and the orbital-phase dependence of the 
interaction, we simulate the interaction between the Be-star 
disc and the neutron star in Be/X-ray binaries, using a 3D 
SPH code. In a general context, such simulations will also 
enable us to study the interaction between the viscous disc 
and the companion in an eccentric orbit. In this paper, which 
is the first of a series of papers dedicated to understanding 
the interaction between the Be disc and the neutron star, we 
study the effects of viscosity on disc truncation in a coplanar 
system. 



2 CALCULATIONS 
2.1 SPH Code 

Simulations presented here were performed with a three- 
dimensional, smoothed particle hydrodynamics (SPH) code. 
The S PH code is based on a version originally developed by 
Benz (Benz 199C; Benz et al. 199C). The smoothing length 



is variable in time and space. The code uses a tree struc- 
ture to calculate the nearest neighbours of particles. The 
SPH equations with the standard cubic spline kernel are 
integrated using a second-order Runge-Kutta-Fehlbe rg in- 



tegrator with individ ual time steps for each particle (Bate 



BonncU & Price 1995), which results in an enormous compu- 



tational saving when a large range of dynamical time-scales 
are involved. 

In our code, the Be disc is modelled by an ensemble of 
gas particles, each of which has a negligible mass chosen to 
be 10"^'' Mq with a variable smoothing length. For simplic- 
ity, the gas particles are assumed to be isothermal at the 
temperature of where T^s is the effective tempera- 

ture of the Be star. On the other hand, the Be star and the 



aeutron star are modelled by two sink particles (Bate et al 
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1995) with corresponding masses. Gas particles which fall 
within a specified accretion radius are accreted by the sink 
particle. We assume that the Be star has the accretion radius 
of 0.9-R, , R* being the radius of the Be star. For the neutron 
star, we adopt a variable accretion radius of 0.9rL, where 
is the Roche-lobe radius for a circular binary. This is because 
having a small accretion radius is time-consuming and it is 
unphysical to adopt an accretion radius smaller than the 
smoothing length of the particles near the disc outer radius. 
An approximate formula for r-L is given by 



: 0.462 



1 + g 



1/3 



D 



(1) 



(e.g., Warner 1995) with the mass ratio q = Mx/M,, where 



Mx and M* are the masses of neutron star and Be star, 
respectively, and the distance between the stars, D. 

The SPH artificial viscosity, Hij , between particles i and 
j enters the momentum equation as 



dt 



^ --^m, + ) V,W{u,,h,,), (2) 

\Pi Pj 



where v is the velocity, m is the mass, P is the pressure, p 
is the density, W is the standard cubic spline kernel, rij is 
the distance between particles i and j, hij — (hi + hj)/2 is 
the mean of the smoothing lengths of particles i and j, and 
Hij has the following standard form. 



{~CesPHCs^^ij + /^SPHMij) /pij '"ij-fij < 



(3) 



( iMonaghan fc Gingold 198^ ), where aspu and /3sph are 
the linear and nonlinear artificial viscosity parameters, re- 
spectively, pij = {pi + pj)/2, Vij = Vi — Vj, Cs is the 
isothermal sound speed, and fiij — hijVij-rij /{rij +rf) with 
rf = O.Ol/i^ 

In the continuum limit, the viscous force in the 
above SPH formalism is written as 



-Fv = 



2p 



[V-{c,pS)+V{c,pV-v)] 



(Meglicki, Wickramasinghe & Bicknell 1993), where 



47r / gdW , 1 
^^^"15 /'^^'^'■^5 



(4) 



(5) 



in the 3D SPH code with the cubic spline kernel and 

is the deformation tensor. If we assume that the density 
varies on a length-scale much longer than the velocity, we 
have 



-Fv = :j^asPHCsft [V^i; + 2V(V ■ v) 



(7) 



This implies that the shear viscosity u and the bulk viscosity 
i^b are given by 



V = — aspHCs/i 
and 



(8) 



(9) 



respectively, in the continuum limit of the 3D SPH formal- 
ism. On the other hand, in the Shakura-Sunyaev viscosity 
prescription, the shear viscosity v is written as 



f = assCsH, 



(10) 



where H is the scale-height of the disc. From equations (|8|) 
and (0), we have the relation between the SPH artificial 
viscosity parameter asPH and the Shakura-Sunyaev viscos- 
ity parameter ass as 

ass = Y^'^^P^Jj^ 

if we assume that \7-v = 0. In general fiows, however, 
V-D 7^ 0. Moreover, the viscosity is artificially turned off for 
divergent flows in our model [see equation (^]. Therefore, 
equation (jllj) should be taken as a rough approximation to 
the relation between ass and aspH. 

In several simulations we report in this paper, we 
adopted osph = 1 and /3sph = 2, in which case ass is 
variable in time and space. In the other simulations, how- 
ever, we adopted constant values of ass in order to roughly 
model the a disc. In these simulations, aspH ~ lOassH/h 
was variable in time and space and /3sph = 0. 

The mass ejection mechanism from the Be star is not 
known. In our simulations, it is modelled by injecting gas 
particles at r = rinj. In our normal resolution simulations, 
which finally have ~ 20, 000 particles, we take Hnj = 1.02ii«, 
while in a high-resolution simulation we performed with 
Qss ~ 0.1, which finally has about 140, 000 particles, we 



take Tin 



1.01-R, . These values were adopted so that 



^inj ~ R* + h/2, where ft is a typical smoothing length near 
the stellar surface. In order to avoid an additional complex- 
ity, we kept the injection rate constant in each simulation. 
Once injected, gas particles interact with each other. As a 
result, most of the injected particles fall onto the Be star by 
losing the angular momentum and a small fraction of par- 
ticles drift outwards, getting the angular momentum from 
the other particles. 

As the Be star, we take a BOV star with Af* = 18Mq , 
Rt = 8Rq, and T^s — 26, 000 K, which has the typical 
spectral type of Be/X-ray binaries. With these parameters, 
the disc scale-height H is 0.02r at r = _R» and increases 
as r^^^. For the neutron star, we take Mx = 1.4Mq and 
Rx = W cm. 



2.2 Initial Configuration 

We set the binary orbit on the x-z plane with the major axis 
along the a;-axis. At t = 0, the neutron star is at apastron. 
It orbits about the Be star primary with the orbital period 
Porb and the eccentricity e. The angle of misalignment, /3, 
between the binary orbital plane and the Be disc plane is 
set to in this paper. The unit of time is Porb, unless noted 
otherwise. 

Each simulation is very time-consuming. It takes three 
to four weeks to perform each of the normal-resolution sim- 
ulations, which runs on a single processor of an HP Exem- 
plar V2500, and about four months to perform the high- 
resolution simulation, which runs efficiently on eight pro- 
cessors of an SGI Origin 3800. This is particularly so for 
long-period systems. Thus, in this paper, we present results 
for only a restricted range of parameter space as the first 
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step in the study of the interaction between the Be disc 
and the neutron star. We fix the orbital period Porb, the 
orbital eccentricity e, and the misalignment angle l3{— 0), 
allowing only the viscosity parameter to vary. We adopt 
Porb = 24.3 d and e = 0.34, targeting 4U 0115+63, one of 
the best studied Be/X-ray binaries. Then, the semi-major 
axis a = 6.6 • 10^^ cm ~ 12i?, and tl given by equation (|l|) 
ranges between 0.13a at periastron and 0.26a at apastron. 



2.3 Testing the Code: Discs around Isolated Be 
Stars 

Before attempting the binary simulations, we applied the 
code to an isolated Be star in order to test whether it gives 
reliable results compared to those obtained by using a one- 
dimensional diff'usion-type equation for the surface density. 



2.3.1 One- Dimensional Model 

As mentioned above, we adopt the viscous decretion disc 
model, which yields a geometrically thin, nearly Keplerian 
disc around the Be star. For simplicity, we assume that the 
disc around an isolated Be star is axisymmetric and Ke- 
plerian and in hydrostatic equilibrium in the vertical {z-) 
direction. The evolution of such a disc is described by 



9E 

dt 

with 

K = - 



ld_ 

r dr 



d ( 2 2.^\ 

or ^ '_ 



(12) 



or ^ ' 



rE 



d_ 

dr 



{r^n) 



(13) 



(e.g., Pringle 198l[ ), where E is the surface density, is 
the radial velocity, and Q — [GM, /r^Y^^ is the angular 
frequency of disc rotation. 

In order to create a situation similar to that in SPH sim- 
ulations, we inject mass at a constant rate at rinj = 1.02_R«. 
At the inner boundary r = Rt, we adopt the torque-free 
boundary condition, E = 0. We also take E = as the outer 
boundary condition at r = 10'^ i?*, which affects the disc 
structure only in a region near r = 10'^ ii*. 

The evolution of E for ass ~ 1 and ass ~ 0.1 for 
the initial three years is shown by dashed lines {t = 1 yr, 
2yr and 3yr from left) in the left panels of Figs. ^ and 
^ respectively. The surface density is measured in units 
of p_iigcm~^, where is the highest local density at 

t = l yr normalised by 10 ~^^gcm~'^, a typical value for Be 



stars (Waters et al. 1988). It should be noted that, in this 
framework, no steady sol ution is prese nt for decretion discs 
unlike for accretion discs (Pringle 1991). Instead, the formal 
solution of equation (^2|) with dT,/dt = and Vr — gives 
the disc structure a,t t ^ oo, which is given by E ~ in 
our isothermal disc model. 



2.3.2 3D SPH Simulations 



Using the model described in §2.1, we performed two 3D 
SPH simulations of the disc evolution around an isolated 
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Figure 2. Evolution of the viscous decretion disc around an iso- 
lated Be star in a 3D SPH simulation with asg = 1. In the left 
panel, the surface density evolution for the first three years is 
shown by the solid lines (t = lyr, 2yr and Syr from left). The 
surface density is measured in units of p-ngcm"'^, where p_ii is 
the highest local density at t = lyr normalised by lO^^^gcm""^, 
a typical value for Be stars. For comparison, the surface density 
distribution at the same epochs in a corresponding ID model is 
shown by the dashed lines. The right panel shows the disc struc- 
ture at t = Syr. The solid, the dashed, the dash-dotted, and the 
dotted lines denote the surface density, the azimuthal velocity 
normalised by the critical velocity of the Be star, the radial Mach 
number, and asPH- In both panels, the density is integrated ver- 
tically and averaged azimuthally, while the velocity components 
and osPH are averaged vertically and azimuthally. The profile 
of V0 is indistinguishable from one proportional to r~^/^. An- 
notated at the bottom of the right panel is the number of SPH 
particles at this epoch. 
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Figure 3. Same as Fig. ti, but for ass = 0.1. 



Be star with qss = 1 and ass ~ 0.1. In these simulations, 
QSPH is variable in time and space and /3sph = to keep 
ass constant, as described in the previous section. 

Figs. H and |^ show the surface density evolution (left) 
and the disc structure at t = 3yr (right) for ass = 1 and 
ass = 0.1, respectively. The time interval between adjacent 
contours in the left panel is 1 yr. In the right panel, the 
solid, the dashed, and the dash-dotted lines denote the sur- 
face density, the azimuthal velocity normalised by the stellar 
critical velocity, and the radial Mach number, respectively. 
In the figures, the density is integrated vertically and aver- 
aged azimuthally, while the velocity components are aver- 
aged vertically and azimuthally. The dotted line shows the 
vertically and azimuthally averaged distribution of qsph re- 
quired to keep ass constant. 

From Figs. |^ and ^, we observe that our 3D SPH code is 
capable of producing the disc evolution similar to that in the 
ID model, except that for small viscosity the disc evolves a 
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0.01 




0.5 ~ 



0.5 



Figure 4. Same as Fig. ^ but for Qsph = 1 and /3sPH = 2. Tlie 
dotted line in tiie riglit panel denotes the distribution of ags, 
after being averaged vertically and azimuthally. 




Figure 5. Evolution of the mass-loss rate from the star, Mdcci 
(thin line) and the disc mass Mjiac (thick line) for a = I (left) and 
a = 0.1 (right). The mass-loss rate is measured at r = R, -f 2rinj 
and averaged over ~ 9 days. 



bit faster than in the ID model. The disc structure is almost 
Keplerian and the radial velocity, which decreases with time, 
is very subsonic for r < lOR, at i > 1 yr. These features are 
in agreement with the observed characteristics of Be stars. 

For the purpose of comparison, we present in Fig. ^ 
the result from a simulation with aspu ~ 1 and /3sph ~ 2. 
In this simulation, ass is variable in time and space and 
is roughly proportional to p~^^^r~^^^ . In the initial phase 
of disc formation, the density distribution has a very steep 
slope in the radial direction, so that ass stays high except 
for a region near the injection radius. As time goes on, ass 
in the outer region decreases, because of the increase in the 
density there. As a result, ass has a local maximum near the 
star, which, in turn, causes a dip in the density distribution, 
the feature not seen in simulations with constant ass- 

Fig. ^ also shows that the disc structure outside the dip 
(r > 3i?*) is, in a rough sense, similar to that in simulations 
with constant ass- Since our purpose is to study the interac- 
tion between the Be disc and the neutron star in a system, 
of which the periastron distance is ~ 8i?*, the difference 
in the disc structure close to the star does not matter very 
much. In the following sections, we will see that both kinds 
of models give similar results. 

In our model, in which the mass injection rate is kept 
constant, the net mass- loss rate from the star through the 
disc decreases as the disc mass increases. This is because 
a larger fraction of injected particles must lose the angular 
momentum and fall back to the star to support a larger disc. 
Fig. ^ gives the change in the mass-loss rate from the star, 
Mdec, (thin line) measured at r = 2rinj — R, and the disc 
mass Md (thick line) for the simulations shown in Figs. ^ 
and 1^. To reduce the fluctuation noise, Mdcc is averaged 
over ~ 9 days. As expected, a larger viscosity makes the 
decrease in A/doc faster. From the figure, we note that, for 
the first several years, Mdcc ~ several x IO^^^P-hA/q yr^^ 
for a wide range of the viscosity parameter. 

It is interesting to compare the model mass-loss rate 
with the observed one. The observed equatorial mass-loss 
rate is, however, a poorly determined quantity. It has been 
measured only during the disc-formation phase for several 
stars. Among them, X Per is the only one star for which both 
the equ atorial mass-loss rat e and the highest disc density are 
known. Telting et al. (1998 ) studied the long-term behaviour 
of the Be disc of X Per (4U 0352-1-30), a Be/X-ray binary 
system with Porb ~ 250 d and e = 0.11, and found that the 
equatorial mass-loss rate for a disc build-up phase of less 



than 230 d was greater than 5.3 x 10 ^Mq yr ^ and that 
for the following 380 d was 3.7 x 10"^AfQyr"\ They also 
found that the base density of the Be disc of X Per was as 
high as 1.5 x 10"^" cm~^. It should be noted that the model 
mass-loss rate shown in Fig. ^ is in good agreement with the 
observed equatorial mass- loss rate of X Per. 

The simulations shown above assume a continuous mass 
supply from the Be star for studying the formation and evo- 
lution of a persistent disc. Some Be stars exhibit a transient 
disc formation by an episodic mass loss. S uch a situation was 
investigated by KroU & Hanuschik (1997). They studied the 
evolution of the gas explosively ejected from a Be star to 
model the transient disc formation and disc decay, using a 
3D SPH code. They found that the gaseous particles form a 
nearly Keplerian disc in the viscous time-scale. Because of 
the difference bet ween the boundary conditi ons adopted in 
this paper and by KroU & Hanuschik (1997), the structure 
of our persisten t disc is different from th at of the decaying 
disc studied by [KroU fc Hanuschik (1997|). 



3 Be DISC-NEUTRON STAR INTERACTION 
IN A COPLANAR SYSTEM 

In the previous section, we have seen that our model with a 
3D SPH code is capable of reproducing the results from ID 
simulations of formation and evolution of the disc around 
isolated Be stars and explaining the observed equatorial 
mass-loss rate from the Be star. In this section, we apply 
our model to a Be/X-ray binary with Porb = 24.3 d and 
e = 0.34, assuming coplanarity between the Be disc and the 
binary. The semi-major axis of the binary, a, is 6.6 x lO^'^ cm 
(~ 12i?, ). We have run simulations with ass ~ 1 for 30Porb 
(~ 2.0 yr), by which time the disc almost reaches an equi- 
librium state, while we have run the other simulations for 
50Porb (~ 3.3 yr). In all simulations but one, the number of 
SPH particles at the end of the simulation was about 2 x 10*. 
In order to confirm the reliability of the results obtained by 
those simulations, we performed a simulation with ass =0.1 
with a higher resolution, in which the number of SPH par- 
ticles at the end of the simulation was about 1.4 x 10^. 



3.1 Disc Evolution under the Influence of the 
Neutron Star 
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Table 1. Summary of the results from binary simulations. The system has Pah = 24.3 d and e = 0.34 in all simulations, 'var' in 
the viscosity-parameter columns means that the quantity is variable in time and space. For the agg = 1 simulation, rd is the radius 
at which the surface density distribution has a break, while for the other simulations it is the radius at which the disc is truncated, 
^dec s-iid Ajfacc are the net mass-loss rate from the Be star and the mass capture rate by the neutron star, respectively, in units of 
P—hMq yr~^, where p_ii is the highest local density at t = lyr normalised by 10~^^gcm~'^, and Si o is the strength of the (1,0) 
mode. < ■ • ■ > denotes the average over t{ — bP^ih £i * £i if • ' — ' in the last column means that no precession is seen. 



Viscosity parameters 
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Artymowicz & Lubow (1994) investigated the tidal/resonant 
truncation of circumstellar and circumbinary discs in eccen- 
tric binaries and found that a gap is always formed between 
the disc and the binary orbit . Follo wing their formulation 



Negueruela fc Okazaki (2001 ) and Okazaki & Negueruela 
(2001a| showed that the Be disc in Be/X-ray binaries is 
truncated via the resonant interaction with the neutron star 
as long as ass <S 1. In the following, we investigate the reso- 
nant interaction between the Be disc and the neutron star in 
more detail, by analysing the results from 3D SPH simula- 
tions. Table |^ lists some characteristic quantities from these 
simulations, which will be discussed below. 

Figs. ^ and ^ show the disc evolution under the influ- 
ence of the neutron star for ass = 1 and ass = 0.1, respec- 
tively. The truncation of the disc is obviously more evident 
for ass = 0.1 than for ass = 1- From Fig. ^, we clearly see 
how the resonant torque works on a viscous decretion disc. 
When the disc size is small {t < lOPorb for ass = 0.1), the 
disc evolution is almost identical to that around an isolated 
Be star. As the disc grows, however, the effect of the res- 
onant torque from the neutron star becomes apparent and 
the radial density distribution begins to have a break at a 
radius around the 4:1 to 5:1 resonance radii (for ass = 0.1). 
We call this radius the truncation radius. Since the reso- 
nant torque prevents disc material from drifting outwards, 
the disc density increases more rapidly than in discs around 
isolated Be stars. Outside the truncation radius, the density 
decreases rapidly. The wavy patterns seen in the surface den- 
sity distribution in the left panel and in the radial velocity 
distribution in the right panel of Fig. |^ are due to the tightly 
wound spiral density wave excited in the disc. 

In contrast to the low-viscosity simulation shown in 
Fig.^, the resonant torque has little effect on the disc struc- 
ture when the viscosity is very high. As seen in the left panel 
of Fig.^, there is an only very modest break in the surface 
density distribution for ass = 1- In this simulation, the disc 
almost reaches an equilibrium state at t ~ 15Porb (~ 1 yr) in 
the sense that the disc structure varies regularly, depending 
on the orbital phase. 

We also performed a simulation with ass = 0.3. The 
disc structure obtained is something between those with 
ass = 0.1 and ass = 1- There is a clear break and a wavy 
feature in the radial surface-density distribution but they 
are not so strong as those for ass ~ 0.1. 

For comparison, we present in Fig. ^ the result from a 
simulation with aspu = 1 and Psph ~ 2, in which ass is 
variable in time and space. From Fig. ^, we observe that the 
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Figure 6. Evolution of the viscous decretion disc with ass = 1 
in a Be/X-ray with Porb = 24.3 d and e = 0.34. The left panel 
shows the surface density evolution. The time interval between 
adjacent contours is 5Porb (~ 1/3 yr) {t = 5Porb, lOPorbi ■■■ 
from left). The right panel shows the disc structure at the end of 
the simulation. The solid, the dashed, and the dash-dotted lines 
denote the surface density, the azimuthal velocity normalised by 
the critical velocity of the Be star, and the radial Mach number. 
In both panels, the density is integrated vertically and averaged 
azimuthally, while the velocity components are averaged vertically 
and azimuthally. The profile of for r < 0.7a is indistinguishable 
from one proportional to r~^^^. Annotated at the bottom of the 
right panel is the number of SPH particles, A^sPHi at this epoch. 
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Figure 7. Same as Fig.H, but for ass = 0.1. 



disc with constant artificial viscosity parameters evolves in 
a similar way to that of the a-disc with a similar viscosity 
parameter, except for the presence of the dip in the den- 
sity distribution close to the star, as was expected from the 
previous section. 
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Figure 8. Same as Fig. Bl but for osph = 1 a^nd PsPH = 2. 
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Figure 9. Same as Fig. ^, but for the high-resolution simulation. 



In order to study the interaction between the Be disc 
and the neutron star in more detail, we performed a high- 
resolution simulation with ass ~ 0.1. In this simulation, 
the number of SPH particles is about an order of magni- 
tude larger than, and so the smoothing length is on average 
less than a half of, that of other simulations. Fig. ^ shows 
the surface density evolution for t — — 45Porb and the 
disc structure at t = 45Porb in this high-resolution simula- 
tion (unfortunately, the allocated computing time ran out 
at t = 47Porb). From Fig. ^, we easily see more detailed disc 
structure than that in the corresponding simulation with 
a lower resolution shown in Fig. M. The wavy pattern in 
the surface density distribution caused by the spiral den- 
sity wave is more remarkable in the high-resolution simu- 
lation. This is because a larger number of particles give a 
higher resolution of the interacting region, which makes the 
interaction more localised and stronger. Although the sur- 
face density profile has breaks near the 5:1 resonance radius 
(r ~ 0.33a) and the 4:1 radius {r ~ 0.39a), the steep den- 
sity decrease already begins at a much smaller radius, which 
coincides with the outermost density peak of the spiral wave. 

In the rest of this paper, we mainly present the results 
from this high-resolution simulation as the representative 
ones with ass ~ 0.1, because it gives a better understanding 
of the star-disc interaction. 



3.2 Phase-Dependent Disc Structure 

Most Be/X-ray binaries with known orbital parameters have 
orbital eccentricities in the range from 0.3 to 0.9. In such sys- 



tems, the star-disc interaction is likely to be strongly phase- 
dependent. In this subsection, we discuss phase-dependent 
features except for the mass capture rate by the neutron 
star, which will be discussed separately. 

Figsjl^ and |l^ give snapshots covering one orbital pe- 
riod for ass = 1 and ass ~ 0.1, respectively. Each panel 
shows the surface density in a range of about two orders of 
magnitude in the logarithmic scale. From these figures, we 
note a remarkable difference in the disc structure between 
ass ~ 1 and ass ~ 0.1. For ass ~ 1, the disc has a signif- 
icant density up to the periastron distance and experiences 
a strong interaction at and after the periastron passage of 
the neutron star. 

On the other hand, for ass = 0.1, the resonant torque 
from the neutron star is much more effective at truncating 
the disc than for the high viscosity disc. The sharp decline 
in the disc density outside the 5:1 resonance causes a gap 
between this radius and the periastron distance, apparently 
reducing the mass capture rate by the neutron star, as will 
be seen in Section 3.4 . 



In both cases, the spiral density waves are clearly seen 
between the periastron passage and the apastron passage. 
The opening angle of the spirals, which is related to the 
effective gravity in the disc, is smaller for a larger binary 
separation. 

Through the resonant interaction, the angular momen- 
tum is transported from the disc to the binary. We have 
to admit that the effect of the angular momentum trans- 
port on the binary turned out to be much stronger than we 
had expected. Despite the fact that the mass of each par- 
ticle is only 10~ ^''Mq so that the disc mass is only about 
10~® — 10~^Mq in our simulations, the increase in the bi- 
nary orbital period is visible in the late stage of simulations. 
The computed phase lags behind the correct orbital phase 
by ~ 0.01 at t ~ 30 in the simulation with ass ~ 1 and 
by ~ 0.04 and ~ 0.05 at t ~ 45 in the high and normal 
resolution simulations with ass ~ 0.1, respectively. In the 
following, figures should be read taking this phase shift into 
account. 

In order to have a measure of the disc radius, i.e., the 
radius at which the disc density has a major break, we ap- 
plied a non-linear least-square fitting method to the radial 
distribution of the azimuthally-averaged surface density E, 
adopting the following simple fitting function. 



E oc 



1 + {r/r^y 



(14) 



where p and q are constants and Vd is the disc radius. 

Fig. |l^ shows the phase dependence of the disc radius 
(thick line) in the simulations shown in Figs.^ and ^ To 
reduce the fluctuation noise, we folded the data on the or- 
bital period over 25 < i < 30 for ass = 1 and 42 < t < 47 
for ass = 0.1. In the figure, the horizontal solid lines denote 
some of the lowest n : 1 resonance radii (the 2:1, the 3:1, . . ., 
the 10:1 from top), while the thin sinusoidal line denotes the 
orbit of the neutron star. The origin of the phase is at the 
periastron passage of the neutron star. 

From Fig. ^ we note that the disc radius coincides 
with the 4:1 resonance radius {r/a — 0.39) for ass = 1, 
whereas the disc has a radius intermediate between the 4:1 
radius and the 5:1 radius {r/a = 0.33) for ass = 0.1 (the 
mean rd is 0.39a for ass = 1 and 0.36a for ass ~ 0.1). The 
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Figure 10. Snapshots for ogg = 1, which cover one orbital period. Each panel shows the logarithm of the surface density. The dark spot 
near the origin is the Be star, while another dark spot orbiting about the Be star denotes the neutron star with the variable accretion 
radius. Annotated at the bottom of each panel are the number of SPH particles, A^spHi and the integrated number of particles captured 
by the neutron star, A^acc ■ 
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Figure 11. Same as Fig.hOl but for the high-resolution simulation with agg = 0.1 



latter is a typical fe ature expected for a disc in which the 
wave damps locally (Artymowicz & Lubow 1994). 



We also note that the disc radius modulates around the 
mean value. The disc shrinks after the periastron passage of 
the neutron star, which gives a negative torque to the disc. 
After that it restores its radius by viscous diffusion, so that 
the amplitude of the modulation is larger for agg — 1 than 
for ass = 0.1. 

Since the disc structure depends on the binary phase, 



the disc emission is expected to exhibit an orbital modula- 
tion as well. To see whether this is the case, we calculated 



the continuum flux from the Be disc for r > 2rin 



suming that the disc is optically thin and the emissivity is 
proportional to p^, where p is the local density. For simplic- 
ity, we ignored the effect of the obscuration of the disc by 
the star, which will become important for systems with high 
inclination angles. We then obtained a base flux curve by 
performing the cubic-spline fitting of the fluxes at apastron 
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Figure 12. Orbital-phase dependence of tlie disc radius r^. Aver- 
aging is done over 25 < t < 30 for ags = 1 (left) and 42 < i < 47 
for ogg = 0.1 (right). The horizontal solid lines denote the 2:1, 
the 3:1, the 4:1, . . ., the 10:1 resonance radii from top to bottom, 
and the horizontal dashed line denotes the phase averaged value 
of rj. The thin sinusoidal line denotes the orbit of the neutron 
star. The periastron passage of the neutron star, which occurs at 
phase 0, is denoted by the vertical dashed line. 



passages. The base flux describes the long-term change in 
the continuum flux. Finally, we computed the orbital mod- 
ulation by subtracting the base flux from the instantaneous 
fluxes. 

Fig. |l^ shows the orbital modulation in the optically 
thin continuum from the Be disc for ass ~ 1 (left) and 
ctss = 0.1 (right). For ass = 0.1, the modulation is more 
clearly seen in the high-resolution simulation (thick line) 
than in the normal one (thin line). In order to reduce the 
fluctuation noise, we folded the data on the orbital period 
over the period annotated in the panel. Contrary to what is 
expected from Fig.^, the continuum exhibits little modu- 
lation for Qss = 1, and about one per cent of modulation is 
seen for ass = 0.1. The negative result for ass = 1 suggests 
that the strongly-perturbed outer disc contributes little to 
the optically thin emission, because of its low density. On the 
other hand, the positive result for ass ~ 0.1, in particular in 
the high-resolution simulation, though the amplitude is still 
very small, seems to come from the region where the den- 
sity is enhanced by the spiral wave, because the disc radius 
does not change signiflcantly and the rise and the subse- 
quent decline of disc emission are in phase with the density 
enhancement and the subsequent damping caused by the 
spiral density wave. The difference between the modulation 
patterns from the high and normal resolution simulations 
with ass =0.1 suggests that resolving the spiral density 
wave is important to obtain a reliable modulation pattern 
of such low amplitude. 

Although Fig. ^ gives little observability of the or- 
bital modulation in optically thin disc emission, it is still 
likely that optically thick emission, such as Balmer lines, 
will show signiflcant orbital modulation. Calculating the op- 
tically thick emission from the disc is, however, beyond the 
scope of this paper. 

3.3 Excitation of the Eccentric Mode 

In a circumbinary disc around an eccentric binary, an ec- 
centric mode is excite d through direct driving due to a one- 
armed bar potential (Artymowicz & Lubow 1996a). In this 
subsection, we study whether the same mechanism works in 
Be/X-ray binaries. We analyse the evolution of the eccen- 
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Figure 13. Orbital modulation in the optically thin continuum 
from the Be disc for agg = 1 (left) and agg = 0.1 (right). In the 
right panel, the thick and thin lines are for the high and normal 
resolution simulations, respectively. Phase corresponds to the 
periastron passage of the neutron star. 



tricity in the Be disc by decomposing the surface density dis- 
tribution into Fourier components which vary as expi{k<j) — 
i^lovht), where k and I are the azimuthal and time-harmonic 
numbers, respectively, and fJorb = [G(A/* + AIx)/a'^]^^^ is 
the frequency of mean binar y rotation. 

Following Lubow (1991), we define the mode strength 

by 



7rMd(l + 5fc,o)(l + 5f,o) 



t + 27rn 



dt' 



X dr dc^^{r,(j>,t)f{k(p)g{lt'), 



(15) 



where / and g are either sin or cos functions. The surface 
density here is computed by summing up 5 functions at par- 
ticle positions, not by taking the kernel into account as has 
been done in the previous sections. Then, the total strength 
of the mode (k, £) is defined by 



'^cos,sin,fc,£ 
~l~'5'sin,cos,A;,^ "t" 'S'sin,sin,fc,^) 



(16) 

Fig. shows the excitation and precession of the (1,0) 
mode (i.e., the eccentric mode) for ass = 1- The upper panel 
shows the strengths of several modes, while the lower panel 
shows the evolution of the angle uJd between the eccentric 
vector of the disc and that of the binary orbit defined by 



tanojd 



(17) 



The eccentric mode grows initially linearly in time 
{t < 8) , as is predicted by the theory and was also found 
by Lubow & Artymowicz (200C). At t ~ 15 the strength of 
the eccentric mode is saturated at Si.o ~ 0.04. As it is sat- 
urated, the mode stops precessing and is locked at Wd ~ i"- 

In order to study the structure of the eccentric mode 
in more detail, we analyse the orbits of individual parti- 
cles. The position and velocity of each particle are instan- 
taneously equal to those of an elliptical Keplerian orbit of 
semi-latus rectum A and the eccentric vector egpH with the 
amplitude egpH and the longitude of periastron , lj, with re- 
spect to the stellar periastron (e.g., Ogilvie 2001 ). The radial 
coordinate r is then given by 
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Figure 14. Excitation of tlie eccentric mode. The upper panel 
shows the strengths of several modes. The solid, the dashed, the 
dash-dotted, and the dotted lines denote the strengths of the (1,0) 
mode, the (2,3) mode, the (1,1) mode, and the (2,2) mode, re- 
spectively. The lower panel shows the angle between the eccentric 
vectors of the disc and the binary. 



A 



1 + espH cos( 



(18) 



Given espH <S 1, A ~ r in our simulations. 

The upper panels of Fig. ^ present the distribution 
of the orbital parameters of disc particles for ass = 1 at 
t = 30. From the left panel, we note that the disc eccentric- 
ity increases almost linearly in A for A < 0.3a and is roughly 
constant for A > 0.3a. The right panel for the radial depen- 
dence of the longitude of periastron shows that the eccentric 
mode is twisted in a trailing sense. 

The normal-resolution simulation with qss ~ 0.1 
showed a similar trend. The eccentric mode grew initially 
linearly in time {t < 40). The mode strength was almost sat- 
urated at t ~ 45 at Si.o ~ 0.1. The mode also stopped 
precessing and was locked at uj^ ~ 3tt/2. 

The distribution of the orbital parameters of disc parti- 
cles at t = 50 from this simulation is presented in the lower 
panels of Fig. |l^. From the left panel, we note that the eccen- 
tricity has a maximum at A ~ 0.3a. The right panel shows 
that the eccentric mode is twisted in a trailing sense, as in 
the case for q ss ~ 1- 

Recently, [Ogilvie (2001 ) has developed a non- linear the- 
ory of the evolution of the shape and surface density of a 
three-dimensional eccentric accretion disc. When the eccen- 
tricity of the disc is small so that terms of relative order 
O(e^) may be neglected, this theory provides a linear evo- 
lutionary equation for the complex eccentricity E{X,t) = 
e(A, f) e'"^'^'*^ of the disc. We have reworked this theory for 
the case of a strictly isothermal decretion disc with no ra- 
dial velocity, and also evaluated the tidal forcing terms as- 
sociated with a companion object of small eccentricity. The 
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Figure 15. Eccentricity espH and the longitude of periastron, 
LU, of the particle orbits from the normal-resolution simulations 
with Qgg = 1 at t = 30 (upper) and agg = 0.1 at t = 50 (lower). 
The gray-scale plot in each panel shows the particle distribution 
in linear scale. < e > and < lu/tt > are mean values of egpH and 
uj/tt. 



governing equation is then of the form 
E(GM*A)^/'^ = ^(ZiEc^^A) + Z2^c^, 



liGMxE/3 
"^4 Yi, 



bl%i(3)E-b]^/^{(3)E^ 



(2) 



(19) 



where E^, and Ab are the (constant) eccentricity and semi- 
latus rectum of the binary orbit, is the Laplace coeffi- 
cient of celestial mechanics, and /3 = A/Ab. The dimension- 
less coefficients Zi and Z2 are given by 



Zi = CiE + C2A-— -, 



Zi^CaE + CiX^, (20) 



where Ci , . . . , C4 are dimensionless complex constants that 
depend only on the shear and bulk viscosity parameters 
ass and ab- In the limit ass,ab these coefficients be- 
come purely imaginary and give rise to precession induced 
by the pressure of the disc. Most importantly, the real part 
of C2, which has the role of a diffusion coefficient for short- 
wavelength eccentric perturbations, turns out to be positive 
when ab/ass > 2/3, as is true in the SPH simulations [see 
equation (^]. If this condition were not satisfied, the disc 
would experience a local eccentric instability and equation 
( p^ could not be evolved forward in time. 

Equation ( p^ ) has the character of a dispersive and dif- 
fusive linear wave equation for E{X,t). The eccentricity of 
the binary provides a tidal forcing that is independent of 
time. Starting from an initially circular state E — 0, the 
eccentricity of the disc first grows linearly in time and then 
approaches a steady state. The steady eccentric shape can 
be determined either by solving the time-dependent problem 
until the transient response decays, or by solving the second- 
order ordinary differential equation obtained by setting the 
time-derivative to zero. We have performed both of these 
calculations, adopting the boundary condition i5 = at the 
stellar surface, a 'stress-free' boundary condition Zi = 
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Figure 16. Model eccentricity e and longitude of periastron uj 
of the particle orbits for ogg = 1 (upper) and ctgg = 0.1 (lower). 
The notional outer edge is at /3(= A/A^,) = 0.39. 



at a notional outer edge (/3 — 0.39), and a surface density 
distribution E oc A~^. 

The steady configurations of e(A) and aj(A) for the cases 
agg — 1 and Qgg =0.1 are illustrated in Fig. y_6[ One should 
not expect a perfect agreement with the SPH simulations 
because the theory is based on the assumption of small ec- 
centricities and only approximates the actual surface density 
distribution. However, the steady eccentric shapes based on 
linear theory offer a fair explanation of the results of the 
SPH simulations, especially for the high-viscosity case. 

Fig. |l^ shows the evolution of an eccentric mode from 
the high-resolution simulation with ass = 0.1. The detailed 
structure of the mode is shown in Fig. ^by snapshots cover- 
ing one precessional period. From these figures, we note that 
the evolution and the structure of the eccentric mode are 
very similar to those from the corresponding normal resolu- 
tion simulation until t ~ 25. At t ^ 25, the eccentric mode 
suddenly begins to precess in a prograde sense, as seen in the 
lower panel of Fig.O. As cOd increases from ~ 37r/2 to 27r, the 
mode strength increases. The mode stagnates at tJd ~ for 
(1 — 2)Porb around t ~ 30, suggesting that the eccentric disc 
at tJd ~ is an unstable configuration. The precession de- 
celerates before reaching lu^ ~ and accelerates after that. 
At the same time, the mode changes its shape. As seen in 
the lower panel of Fig. |l^, the eccentric mode is twisted in 
a trailing sense for ti;d < and in a leading sense after that. 
At ujd ~ 0, apsidal axes of all particle orbits align. 

As seen in Fig. |l^, the twist of the mode increases, as the 
apsidal axis of the disc moves toward the stellar periastron 
(cod — n). The strong radial dependence of the phase 
caused by the twist reduces the eccentricity of the disc, as 
shown in the upper panel of Fig. ^ The precession is fastest 
at u>d ~ TT, suggesting that there is a stable point at ujd ~ 
TT. After that, the mode becomes trailing and its strength 
increases. The precessional period is about 20Porb- 

It is important to note that a similar behaviour is found 
in circumbinary discs around ec centric binaries. According 
to Lubow fc Artymowicz (2000| ), this behaviour occurs as 




Figure 17. Same as Fig. |l^ but for the high-resolution simula- 
tion with ogg = 0.1. 



follows; When the eccentricity of the disc edge is small, ujd 
is locked at a stable value cud = Sn/2. However, as the eccen- 
tricity grows, the locking action weakens, and the prograde 
precession due to the quadrupole moment of the binary po- 
tential dominates. The disc edge begins to precess when its 
eccentricity becomes (0.2 — 0.7)e, and afterwards the eccen- 
tricity oscillates with a precessional period. The disc typi- 
cally attains the eccentricity of (0.5 — l)e. 

We note that the above behaviour of the eccentric mode 
in circumbinary discs around eccentric binaries is strikingly 
similar to that shown in Fig. |l^ except that, in our sim- 
ulation for Be/X-ray binaries, the growth and precessional 
time-scales of the eccentric mode are much shorter and the 
disc eccentricity attained is significantly smaller than in cir- 
cumbinary discs. 



3.4 Mass Capture Rate by the Neutron Star 

As mentioned in Section |l|, most of the Be/X-ray binaries 
show transient X-ray activities. Among them, some exhibit 
periodical X-ray outbursts called Type I, which coincide 
with the periastron passage, while the others show occa- 
sional giant X-ray outbursts called Type II and little or no 
detectable X-ray emission in quiescent phase. 4U 0115-1-63, 
the system we are modelling in this paper, belongs to the 
latter group. In this subsection, we first study how much 
mass is captured by the neutron star in a general context, 
and then discuss whether our model predicts the mass cap- 
ture rate consistent with the observed X-ray behaviour of 
4U0115-f63. 

Fig. |l9| shows the change in the mass capture rate by 
the neutron star, Mace, and the disc mass Md for ass = 
1. The upper panel shows the evolution of Mace and Md, 
while the lower panel shows the orbital-phase dependence 
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Figure 18. Change in the eccentricity egpjj and the longitude lu of the disc particle orbits in the high-resolution simulation with 
ogs = 0.1. Epochs are chosen for illustrative purpose, so the time interval is not constant. The gray-scale plot in each panel shows the 
particle distribution in linear scale. < e > and < uj/it > are mean values of egpH S'lid t^/f at each epoch, respectively. 



of Mace- In the lower panel, we folded Afacc on the orbital 
period over 25 < t < 30 to reduce the fluctuation noise. 
The horizontal dashed line and dash-dotted line denote the 
mean mass capture rate by the neutron star and the mean 
mass-loss rate from the Be star, respectively. 

We have already seen that there is little truncation 
of the Be disc for ass = 1- Fig- confirms this. For 
t > 20 the Be disc is almost in equilibrium in the sense 
that the disc mass and the mass capture rate only shows 
regular orbital modulations with constant amplitude. For 
25 < t < 30, the mean mass- loss rate from the Be star is 
2.4 X 10~^''p_iiMQyr~^ , while the mean mass capture rate 
for the same period is 2.1 x lO^^^p-nMQyr"^. Thus, in this 
high viscosity disc, the neutron star captures the disc mass 
at about the same rate as the mass-loss rate from the Be 
star. 

The simulation with ass = 0.3 also gives a high fraction 
(about 60 per cent for 45 < t < 50) of mass capture rate 
with respect to the mass-loss rate, as listed in Table |l[ This 
indicates that the resonant truncation effect is not effective 
for ass = 0.3, either. 

In contrast to the simulations with ass = 1 and ass ~ 
0.3, the high-resolution simulation with ass = 0.1 revealed a 
subtle feature in the mass capture rate, as shown in Fig. |2o[ 
Before the precession of the eccentric mode began, the mass 
capture rate Mace increased monotonically: Mace = before 
t 10 with the resolution of this simulation. Then, the 
peak mass capture rate showed a gradual increase to a level 
at Mace = (8 - 9) X lO""p_iiM0 yr"^ at t ~ 25. 

After the eccentric mode began to precess. Mace exhib- 
ited a long-term modulation, depending on the longitude 
of the eccentric mode, Ud- It gradually increased as Ud in- 
creased from to 7r and decreased as Ud increased from tt 
to 27r. At the periastron passage at t ~ 38, Mace was a max- 
nnum of (2 - 3) X lO~^V-iiA^0yr"\ which was about a 
factor of three higher than the level before precession. The 
eccentricity of the disc was ~ 0.1 at Wd ~ and 0.02 — 0.03 
at ujd ~ TT. It should be noted that even this small eccen- 
tricity caused a factor of three enhancement in the mass 
capture rate. If the disc had a much stronger disturbance, 
the enhancement in Mace is likely to be much larger. 
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Figure 19. Evolution of the disc mass and the mass capture rate 
by the neutron star (upper) and the orbital-phase dependence of 
the mass capture rate (lower). In the upper panel, the thin line 
denotes the mass capture rate Mace and the thick line denotes 
the disc mass Mj. In the lower panel, the data are folded on the 
orbital period over 25 < t < 30. The horizontal dashed line and 
the dash-dotted line in the lower panel denote the mass capture 
rate by the neutron star and the mass-loss rate from the Be star 
averaged over the period annotated at the bottom of the panel. 

In addition to the long-term modulation due to the pre- 
cession of the eccentric mode, we note that the mass cap- 
ture rate by the neutron star is much smaller and more 
strongly phase-dependent for ass = 0.1 than for ass ~ 1- 
For ass = 0.1, the phase at which the neutron star cap- 
tures the disc mass most slightly lags behind the perias- 
tron passage, because of the presence of a gap between the 
disc outer radius and the periastron. In the high resolution 
simulation, the peak mass capture rate for 42 < f < 47 
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Figure 20. Same as Fig. |l9| but for the high-resolution simula- 
tion with ass = 0.1. In the lower panel, the data are folded on the 
orbital period over 42 < t < 47. For comparison, the mass cap- 
ture rate from the normal-resolution simulation with agg = 0.1 
(thin solid line) is overlapped in the lower panel. 



is ~ 10~^''p_iiAfQyr~^, which is one order of magnitude 
smaller than that for ass = 1- The mass capture rate then 
decreases by two orders of magnitude by the apastron pas- 
sage. Even in the normal resolution simulation, in which 
the disc density around the truncation radius is significantly 
higher than that in the high-resolution simulation, the peak 
mass capture rate is about a factor of three smaller than 
that for agg — 1, and decreases by a factor of fifty by apas- 
tron. Note that similar, stro ngly phase- dependent accretion 
was found in simulations by Artymowicz & Lubow (1996h) 



for circumbinary discs around eccentric binaries. 

In the high-resolution simulation with agg = 0.1, the 
neutron star, on average, captures the disc mass at the 
rate of 2.3 x lO""p-iiA/0yr"^ for 42 < t < 47, while the 
mean mass-loss rate from the Be star for the same period 
is 1.6 X lO~^''p_iiM0yr~^ . This indicates that, even after 
three years of disc growth, about 6/7 of the gas lost from 
the star still accumulates in the disc, making the disc con- 
tinually denser. 

The transient nature in the X-ray activity of Be/X-ray 
binaries is considered to result from the interactions between 
the a ccreted material an d the rotating magnetised neutron 
star (Stella et al. 1986). If the accreted material is dense 
enough to make the magnetospheric radius smaller than the 
corotation radius, the accretion onto the neutron star causes 
a bright X-ray emission (the direct accretion regime) . Other- 
wise, the magnetospheric radius is larger than the corotation 
radius, and the accretion onto the neutron star is prevented 
by the propeller mechanism. The system is then in quies- 
cence (the pr opeller regime). 

Recently, [Campana et al. (200l| ) found that 4U 0115-f 63 
showed a dramatic X-ray luminosity variation from Lx 

X 



close to periastron, whereas it showed a nearly constant X- 
ray luminosity at a level of Lx 2 x lO'^'^ergs"^ near apas- 
tron. They concluded that the system was in the transition 
regime between the direct accretion regime and the propeller 
regime close to periastron and in the propeller regime near 
apastron, because the direct accretion regime and the pro- 
peller regime, respectively, apply for Lx ^ C^'^'^lO^'^erg s~^ 

and L^<^''/'^2 x 10^''ergs"\ where 0.5 < ^ < 1 is a param- 
eter which depends on the geometry and physics of accre- 
tion. 

Below we try to compare our simulat ion results with 
the above criteria by Campana et al. (2001). Although it is 



2 X 10^*ergs-i to L 



5 X 10^'' ergs ^ in less than 15 hr 



likely that the captured material will form an accretion disc 
around the neutron star, we do not know how the accretion 
rate is related to the mass capture rate shown in Figs. |l^ and 
p^ , which is based on the moments at which particles enter 
the variable accretion radius of the neutron star. Therefore, 
we consider two extreme situations, in which tacc ^ Porh 
or tacc ~ L'orb, where tacc is the accretion time-scale. In the 
former situation, the accretion rate profile is approximately 
the same as the profile of mass capture rate, whereas in the 
latter situation, the variation in the accretion rate will be 
much smaller than that in the mass capture rate. We assume 
that the X-ray luminosity is given by Lx = rjGMx M / Rx 
with rj — 1. 

For agg = 1, Lx ~ IQ^'^P-hMq yr~^ at periastron if 
^acc <C L'orb and all the captured mass accretes onto the neu- 
tron star. Note that this level of luminosity enters the direct 
accretion regime. Since 4U 0115-1-63 is considered to be in 
the direct accretion regime only in occasional giant X-ray 
outbursts, we conclude that the agg = 1 model is inconsis- 
tent with the observation if tacc <C L'orb. If iacc ~ L'orb, the 
neutron star will emit the X-ray luminosity corresponding 
to the mean mass capture rate of ~ 2 x lO~^°p_iiM0 yr~^ 
(see the bottom panel of Fig. p^ ). The X-ray luminosity is 
then Lx ~ 2 x lO'^^p^nerg s"'^. This level of luminosity en- 
ters the transition regime. Hence, the ass ~ 1 model is not 
ruled out by the observed constraint if tacc ~ L'orb. 

We have the same conclusion for the agg — 0.3 model. 
Since it gives the mass capture rate about a half of that for 
Qss ~ 1, the model is inconsistent with the observation if 

tacc < Lorb, but is UOt rulcd out if tacc ~ L'orb. 

On the other hand, for agg — 0.1, Lx ~ 
IO'^^P-hMq yr~^ at periastron even if tacc ^ L'orb and all 
the captured material accretes onto the neutron star. This 
suggests that the system is in the transition regime even 
at periastron. It is likely that the low mass capture rate in 
this simulation put the system into the propeller regime in 
most of the orbital phases. Thus, the agg = 0.1 model for 
4U 0115+63 is consistent with the observed X-ray behaviour, 
irrespective of the accretion time-scale. 



4 SUMMARY AND DISCUSSION 

In this paper, we have presented results from 3D SPH sim- 
ulations of the disc formation around isolated Be stars and 
of the interaction between the Be star disc and the neutron 
star in Be/X-ray binarie s, based on the viscous decretion 
disc model for Be stars (Lee et al. 1991). In several sim- 
ulations we adopted constant values of artificial viscosity 
parameters agpH and /3gpH, for which the Shakura-Sunyaev 
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viscosity parameter ass is variable in time and space. In 
the other simulations, we adopted constant values of ass, 
for which aspH is variable in time and space and /Jsph = 0. 
We preferred constant ass simulations because the analysis 
of the results becomes easier. In all simulations, the Be disc 
was nearly Keplerian and the radial velocity was highly sub- 
sonic. These features are consistent with those observed for 
Be discs. 

The simulations of isolated Be stars showed that our 
code is capable of producing results similar to those from 
the ID simulations. The simulated mass-loss rate from the 
Be star in the first several years of disc formation was several 
X 10~^°P_iiMq yr~^ for a wide range of viscosity parame- 
ter, which is consistent with the observed equatorial mass- 
loss rate. Here, p_ii is the highest local density at f = lyr 
normalised by 10~^^gcm~'^, a typical value for Be stars. 

In binary simulations, we have studied the effect of vis- 
cosity on the star-disc interaction in the case of a coplanar 
system with Porb = 24.3 d and e = 0.34, the parameters for 
4U 0115-1-63, one of the best studied Be/X-ray binaries. We 
have chosen these orbital parameters because they enable 
us to compare the simulation results with the observed ones 
and the short orbital period makes the computing time bear- 
able. Some of the results from these simulations are sum- 
marised in Table 0. 

Our simulations showed that there is a radius outside 
which the azimuthally averaged surface density decreases 
steeply. For a smaller ass, the slope outside this radius 
is steeper, giving a stronger truncation effect on the disc. 
Among the simulations with ass = 1, 0.3, and 0.1, we found 
that the resonant truncation of the Be disc is effective only 
for ass = 0.1. For ass = 1 and 0.3, the neutron star cap- 
tures the disc mass at a rate comparable to the mass-loss 
rate from the Be star. These results confirm the pre vious 
semi-analytical result by Negueruela fc Okazaki (2001 ) and 



Okazaki & Negueruela (2001a) on disc truncation that the 



resonant truncation is effective for a disc with ass 1. The 
truncation radius for ass =0.1 roughly agrees with that 
derived semi-analytically. 

Our simulations, in particular the high-resolution simu- 
lation with ass = 0.1, showed how the disc grows under the 
influence of the neutron star. In the initial build-up phase, 
the disc evolution is similar to that for isolated Be stars. But, 
later on, the effect of the resonant torque becomes apparent, 
preventing the disc gas from drifting outwards at several res- 
onance radii. The effect is most remarkable at the 4:1 and 5:1 
radii (r/o ~ 0.39 and 0.33, respectively) for ass = 0.1. As a 
result, the disc density increases more rapidly tha n that for 
isolated Be stars. This feature is consistent with Zamanov 
et al. (^001), who found that the discs in Be/X-ray binaries 



are about twice as dense as those of isolated Be stars. 

Since the neutron star orbits in an eccentric orbit about 
the Be star, the interaction is phase dependent. The disc 
shrinks a little at periastron and then recovers gradually. 
Consequently, the disc emission will vary with phase. For 
our system with e = 0.34, this effect turned out not to be 
large enough to cause an observable variation in the disc 
continuum as long as it is optically thin. It is possible, how- 
ever, that the orbital modulation in the disc continuum in 
systems with much higher orbital eccentricity is observable. 
Moreover, it is likely that the profiles of optically thick emis- 
sion lines from the disc, such as Balmer lines, will show an 



orbital modulation even in systems with moderate eccentric- 
ity, because the disc structure is made non-axisymmetric by 
the periodic excitation and damping of the spiral density 
wave. 

In the Be disc in Be/X-ray binaries, an eccentric mode 
is excited through direct driving due to the (1,0) harmonic of 
the binary potential. In high-viscosity simulations, the mode 
grows initially linearly in time and then approaches a steady 
state, of which steady eccentric shapes based on linear the- 
ory of the evolution of a 3D eccentric decretion disc de- 
scribed in Section 3.3 offer a fair explanation. The strength 



of the mode in the steady state is larger for a smaller value 
of ass. 

On the other hand, in the high-resolution simulation 
with ass = 0.1, the eccentric mode undergoes prograde 
precession at some point. The precession period is about 
20Porb- The precession rate is not constant. It accelerates 
for < ciJd < 71" and decelerates for tt < tJd < 27r, where ujg^ 
is the angle between the eccentric vector of the disc and that 
of the binary orbit. The precession rate is radius-dependent. 
It is larger at a larger radius. Since the eccentric mode is 
leading for < Wd < ti" and trailing for tt < tJd < 27r, the 
twist of the mode increases with time for < ojd < tt and 
decreases for -n < uid < 2ti. The strength of the eccentric 
mode, which anti-correlates with the twist of the mode, de- 
creases for Q < uid < IT and increases for tt < ojd < 27r. 

As the mode precesses, the mass capture rate, A/acc, by 
the neutron star modulates. It is a maximum at ~ tt. 
Even for an eccentricity of 0.1, Mace at maximum is about 
a factor of three higher than the level before precession. If 
the disc had a much stronger disturbance, the enhancement 
in Mace is likely to be much larger. Such a system may tem- 
porarily show periodic X-ray outbursts. 

In our model for Be/X-ray binaries, in which the disc 
material overflows toward the neutron star via the L\ point, 
the mass capture rate by the neutron star becomes strongly 
phase dependent. The dependence is stronger for a smaller 
value of viscosity. In the disc with ass = 0.1, the mass cap- 
ture rate decreases by two orders of magnitude between pe- 
riastron and apastron passages, and after apastron passage, 
no disc mass is captured by the neutron star. Note that our 
model gives a much stronger contrast in the mass capture 
rate than that expected for the stellar wind accretion. 

We have also compared the simulated mass capture 
rate with the observed X-ray behaviour of 4U 0115-)-63, con- 
sidering two extreme situations, in which face ^ Porb or 
tacc ~ Porb, where tacc IS the accretion time-scale. We found 
that the disc model for ass =0.1 gives a result consistent 
with the observation for both situations, whereas the higher 
viscosity models for ass = 1 and 0.3 are ruled out unless 

^acc ^ Porb • 

Analysing multi-wavel ength long-term monitor ing ob- 
servations of 4U 01154-63, Negueruela et al. (2001) found 
that the Be star undergoes quasi-cyclic (~ 3 — 5yr) activ- 
ity, losing and reforming its circumstellar disc. They also 
found that, at some point, the growing disc becomes unsta- 
ble to warping and then tilts and starts precessing. Type 
II X-ray outbursts take place after the warping episode. As 
shown in this paper, our ass = 0.1 model explains many 
of the observed features of 4U 0115-1-63 in the phase before 
the warping occurs. Our simulation, however, showed no dy- 
namical instability, and therefore is incapable of explaining 
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the warping episode. Including the effect of radiation from 
the Be star may turn out to be essential to have a model 
that explains the whole cycle of the disc evolution, which is 
beyond the scope of this paper. 

In this paper, wo have concentrated our study on the 
viscous effects on star-disc interaction in a coplanar system 
with fixed orbital period and eccentricity. This was done not 
only as a first step to have a comprehensive understanding 
of the interaction between the Be disc and the neutron star 
in Be/X-ray binaries, but also to have an archetypal model 
with which we can compare the results from our future sim- 
ulations. In the next paper, we will study the effects of the 
misalignment angles between the Be disc and the binary. 
We will discuss the effects of the orbital eccentricity in the 
third paper, and the fourth paper will conclude the series 
by studying the effects of the orbital period. 
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